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We present a general quantum circuit design for finding eigenvalues of non-unitary matrices on 
quantum computers using the iterative phase estimation algorithm. In particular, we show how the 
method can be used for the simulation of resonance states for quantum systems. 



I. INTRODUCTION 

In the circuit model of quantum computation, a con- 
trolled quantum system, considered to be a quantum 
computer, propagates from one state to another through 
the application of a local propagator called a quantum 
gate. A given algorithm, or computation, can be im- 
plemented on a quantum computer using a sequence of 
quantum gates. These gates can be represented by uni- 
tary matrices in the computational basis. While this 
unitary representation is adequate for many problems, 
it impedes application of quantum computing to new 
problems, where computational bases can only be de- 
scribed through non-unitary matrices. Furthermore, any 
attempt at approximation of these non-unitary parts neg- 
atively impacts the accuracy of the results. For instance, 
the algorithm of Harrow et. al |T] for solving linear sys- 
tems of equations requires a non-unitary quantum gate. 
The circuit designed for this algorithm [2] uses an ap- 
proximated gate for the non-unitary parts. 

Quantum algorithms are known to be more efficient 
than their classical counterparts. One of these algorithms 
is the phase estimation algorithm [3H5], used for finding 
the eigenvalues of a given unitary matrix, or equivalently, 
the eigen-energies of the corresponding quantum system. 
The algorithm has been applied in quantum chemistry 
to obtain energies of molecular systems [6l-[l0] . The algo- 
rithm has also been demonstrated experimentallv [TTl 1 1 2 j . 
However, the simulation of non-hermitian matrices in 
quantum theory requires the simulation of non-unitary 
evolution matrices. 

It has been shown that a quantum circuit can be gen- 
eralized to a non- unitary circuit, whose constituents arc 
non- unitary gates representing quantum measurement. 
Furthermore, it is shown that a specific type of one-qubit 
non-unitary gates, the controUed-NOT gate, and all one- 
qubit unitary gates constitute a universal set of gates for 
the non-unitary quantum circuit without the necessity 
of introducing ancilla qubits [TS]. Recently, Nori et al. 
[14j . proposed a measurement-based quantum algorithm 
for finding eigenvalues of non-unitary matrices. Their 
method draws on ideas from conventional phase estima- 
tion algorithm, frequent measurement, and techniques in 
one-qubit state tomography. 

In this paper, we introduce a systematic way of esti- 



ancilla 



|0) — H 



|0) — H 



mam 



|0) >^ 



H — 



H — 



FIG. 1: Universal circuit design illustrated for two 
qubit operators: The uniformly controlled gates in the 

middle form all elements of U on the diagonal of V. 
The Hadamards at the end carry the same row elements 
to the first rows of each Vi. The initial Hadamards and 
the SWAPs modify the input. 



mating the complex eigenvalues of a general matrix us- 
ing the iterative phase estimation algorithm with a pro- 
grammable circuit design [15) . Our model can also be 
used as a general circuit equivalent for non-unitary ma- 
trices. We also present the application of our method to 
the example of non-hermitian Hamiltonians, which ex- 
hibit predissociating resonances. Since the circuit design 
is universal, it can be used to estimate the eigenvalues of 
any given system. 



II. UNIVERSAL CIRCUIT FOR NON-UNITARY 
MATRICES 



Recently |15j . we described a universal programmable 
circuit design that can be used to simulate any given 
matrix by simply determining the angle values for gates 
using matrix elements. This universal circuit, illustrated 
in Fig{T]for two qubit systems, is constructed as follows: 
For an arbitrary n-qubit system represented by an iV x 
matrix U (where A^ = 2"), (n-l- 1) ancilla qubits are used 
to control N x N block quantum operations in order to 
generate circuit equivalences for rows of U independently. 
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as the output we get from J7|f/') = with a normal- 
ization constant k (the dashes are used to separate the 
anciUa and the main qubits). Next, the construction of 
V and Yp) in Eq.Q are explained in following three steps: 
Formation, Combination, and Input Modification. These 
three steps result in the circuit shown in Fig[T] 



FIG. 2: Each block operations forms a different row of 
the given matrix as their leading row. 



Formation step: 



This is shown in Figj2] 

The matrix representation of Figj2] V, has N block 
matrices located on its diagonal. Each of these blocks, 
Vi , generates the ith row of the given matrix as their first 
rows (the elements of V outside of the blocks are zero) . 

The relationship between an arbitrary input, \Tp), to 
the system U, and the corresponding output, |/3), can be 
defined as U = \(3). Using a modified input, \ip), the 
system with the ancilla, V, can also be used to simulate 
the application of U on predetermined N output states 
obtained from V\(p): 
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in which the states: k^i|0 ... 000-0 ... 0), k/3i|0 . . . 010- 
. . . 0), . . . , kPn]! . . . 110 — ... 0), are exactly the same 



Each element of U is formed using one rotation gate 
Rj located on the diagonal of Vf. Vf corresponds to a 
2iV^ X 2iV^ matrix representing the Formation step: 



Vf 



cos{^) sin{-^) 
—sin{^) cos{^) 



(2) 

where Rj creates the j'th-element, Uj, of U by the equal- 
ity cos(^) = Uj. Here, we have used rotations around 
the y-axis. However, for the complex elements of J7, 
a combination of R^ and Ry gates are needed for each 
Rj. The circuit equivalent of the matrix Vf constitutes 
N"^ uniformly controlled rotation gates (the total number 
of elements of U) operating on the last qubit ((2n+l)st 
qubit). 



2. Combination step: 

In this step, we bring the same row elements of U to 
the leading rows of each Vi by applying a combination 
matrix Vc to Vf, formed in the previous step: V = VcVf. 
This is described as: 
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where each Vc has the following form: 
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For a system of {n + 1) qubits, the matrix in Eq.Q rep- 
resents the Hadamard gates on the first n qubits with 
k = l/-\/2". Thus, for a system of 2n -I- 1 qubits, Vc can 
be defined as the application of the Hadamard gates to 



(n + l)th, nth, 
in the circuit. 



3rd, and 2nd qubits from the bottom 



3. Input modification step: 

After forming all the rows of U through the leading 
rows of each block Vi, we transform the input so that the 
resulting output includes the simulation results defined in 
Eq. ([T]) . The initial input to the final matrix V in Eqjsjcan 
be defined as lip) = |0..0) lijj) where |0...0) is the input 
to the ancilla qubits. In the application of V, we also 



2 



phase 



mam 



m- H 



R-Awk) - H - 



FIG. 3: The iterative phase estimation algorithm for 
the fcth iteration [^-(5]. In the circuit jV'j) is an 
eigenvalue of U, and the angle Wk of the gate 
depends on the previously measured bits defined as 
Wk = -2Tr{0.0xkXk-i.--Xm)binary, where m is the 
number of digits determining the accuracy of the phase 
(pj . Note that Wk is zero in the first iteration of the 
algorithm. 



need to disregard the effects of the elements represented 
by "*" between kuij and kui(^j+iy If the corresponding 

elements in the input are made zero, as in Eq. 
these elements are not-effective in the application of V 
to the following input state: 



\ip) =[ai a2 . . . QfAT . . . 0]^ 
[nai Ka2 ... KCkjv ■ • 
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where k is a normalization constant. This modification 
can be achieved by applying Hadamard gates on the first 
n qubits, and sequential swap operations between the 
(n + l)th and the remaining last n qubits. 

Consequently, the combination of these three steps 
yields the final circuit, which is shown in Figjl] for a 
two qubit system. This circuit can be used to simu- 
late any matrix U having matrix elements less than or 
equal to 1. Finally, the following normalized set of N 
states exactly simulates the given unitary U: |0 . . . 000 — 
... 0) , ... 1 1 ... 110 — ... 0) . For simplicity, in the phase 
estimation algorithm, we will swap the first n ancilla 
qubits with the main qubits and use the following states 
instead: 
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III. FINDING EIGENVALUES OF 
NON-UNITARY MATRICES 

The polar form of a complex eigenvalue, Xj, belonging 
to a matrix, U, can be written as: 



A, 



(7) 



where (j)^ is the phase and = \Xj\ is the magnitude 
of the jth eigenvalue. For unitary matrices, since |Aj| 



is equal to 1, the phase estimation algorithm (PEA) is 
a well-known efficient quantum algorithm used to esti- 
mate the phase, 4>j, and consequently, the corresponding 
eigenvalue, Xj. The iterative phase estimation algorithm 
IPEA [Mi] shown in Fig[6] can also be used to estimate 
the value of the phase, (pj. In each iteration of the phase 
estimation process, we estimate one bit of the m-digit 
binary expansion of the phase represented as: 
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where m is the total number of iterations, which deter- 
mines the precision, x/c is the binary value found from 
the fcth iteration of the IPEA using the 2™~'^th power of 
the matrix, U"^ 

However, for non-unitary matrices, the phase estima- 
tion algorithm alone is not enough. This is because two 
parameters, 0j and rj, must be found to compute the 
eigenvalue. In the following sections, we show that by 
applying the iterative version of the phase estimation al- 
gorithm to the general circuit described in the previous 
section, parameters (pj and Vj can be determined, and 
hence eigenvalues of non-unitary matrices can also be 
found by the phase estimation algorithm without addi- 
tional cost. Since our circuit design has a fixed size and 
scheme, the application of IPEA to the circuit design 
also has a fixed design. Therefore, we only need to set 
the angle values in each iteration of IPEA. 



A. Estimation of <pj 

To use the general circuit design shown in Figjl] for 
4x4 matrices, within the IPEA framework, we need 
to determine the relationship between the chosen states 
and the phase qubit resulting in the measurement. As 
described in the previous section, the circuit scales the 
expected output of the chosen states by k, which clearly 
originates from the Hadamard gates. In phase estima- 
tion, since all of the gates are controlled by the upper 
qubit, this scaling effect of k exists only when the phase 
qubit is one. In order to have the same scaling when the 
phase qubit is zero, we allow the Hadamard gates in the 
circuit to operate without any control qubit. Therefore, 
controlling the swap gates at the beginning and at the 
end of the circuit do not affect the chosen states. Hence, 
they are also allowed to operate without the control on 
the phase qubit. The resulting circuit is shown for two 
qubit systems in Fig|4| Here, we control all operations 
except the Hadamard and the swap gates by the phase 
qubit. 

We estimate the phase on the selected N states, where 
the ancilla is zero in the output. In one iteration of 
the phase estimation algorithm, the effect of the eigen- 
value, represented in polar form, can be observed in 
the evolution of the system as follows: After the first 
Hadamard gate on the phase qubit, as seen in the circuit 



3 



in Fig.(|6|, and the controlled operation (as stated above, 
the Hadamard and the swap gates are not controlled), 
we have the following wave function describing the cho- 
sen states in the first iteration (the ancilla is ignored since 
it is zero in the chosen states): 

vI/i = -^(|0)|^,)+e^2'^(°-^)|A,||l)|^,)), (9) 

where \'ipj) is the input, an estimate of the eigenvector, 
and the first qubit is the phase qubit. Application of 
the rotation gate around the z-axis on the phase qubit 
in the circuit guarantees that only one bit, in Q.Xk, of 
the phase is estimated in each iteration. After the last 
Hadamard applied to the phase qubit, the state in Eq.([9]) 
becomes: 

- ^((^)l^.> +e-(°-)|A,|(MHl)),^^.^) 

(10) 

or more concisely: 

*2 = ^((l + e"-(°-^)|A,|)|0) + (l-e*2'^(°-^)|A,|)|l))|^,). 

(11) 

When xi — 0, since e'-'^'^^ = 1, the above equation reduces 
to: 



2= 2((1 + 1^.1)10) + (1-|A,I)|1))|V',) 



In the case when xi = 1, Eq.(ll) becomes: 



Since f (1 - \Xj\) < f (1 + \Xj\); if the phase bit is 0, 
then xi = 0. If the phase bit is 1, then the xi = 1. 
Thus, the measurement on the phase qubit determines 
xi . In each iteration of the phase estimation algorithm, 
a different bit of the phase is determined in this manner. 
For the fcth iteration of the phase estimation algorithm 
to determine the fcth digit of the phase, we see f (1 — 

|Ajf ') < 1(1 + |Ajf '). Hence, for |Aj| < 1, the more 
we increase the power of U or the accuracy, the more the 
amplitudes of the states, |0) and |1), come closer. 



B. Computing |Aj| 

As we have shown above, the output probability of the 
phase qubit when the ancilla is zero is determined by 
f ((1 + \Xj\)- We measure |0) or |1) with probability P 
equal to |((1 + |Aj|). Hence, we get 
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FIG. 4: The iterative phase estimation algorithm with 
the circuit shown for a 2 x 2 matrix. The decomposition 
of the middle multi control network requires 16 CNOTs 
and 16 single rotation gates. 



can be determined from Eq.( 14). Therefore, the accuracy 



of the estimate value of |Aj| can be further improved by 
using all iterations of the phase estimation algorithm: 
For the /cth iteration, as a general form, the following is 
obtained: 



P, 



W = (f((l + |A,f")f 
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- 1. 



Taking the average of the estimates from different phase 
estimation iterations, the estimate of |Aj| may become 
(12) more accurate. Finally, by using values of |Aj| and (pj 
in Eq. ([t]), we can compute the eigenvalue of the non- 
unitary matrix. 

Note that the complexity of computing the power of 
matrix U on classical computers is O(mM), where m is 
the number of iterations and M is the complexity of ma- 
trix multiplication (this is estimated using the successive 
squaring method). 



IV. APPLICATION TO NON-HERMITIAN 
QUANTUM SYSTEMS 



Suppose we have the non-unitary operator U = 
^—Mt/h ^ non-hermitian Hamiltonian H with energy 
eigenstates \4>j) and corresponding energy eigenvalues 
Ej, i.e., Hlipj) = Ejltpj). Since Ej is an eigenvalue of 
^; if i and h are set to 1, then e"*^^ is the eigenvalue 
of U. Therefore, a N x N non-unitary transformation U 
has eigenvectors \(pi), IV2), I'Pn), with corres pon ding 



eigenvalues A^ 



As shown in Sec 
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eigenvalue \j of the matrix U can be estimated 
also allows us to determine the corresponding eigenvalue 

e 



Ej of the Hamiltonian Ti. from Aj 



P = 

I A, I 



(2((1 + |A.I))', 

2^/P 



(14) 



Since k and P are known (for a 2 x 2 matrix, the value 
of K is 2 due to two Hadamard gates in the circuit), |Aj| 



A. Example: Simulation of Resonance States 

As an example, we consider the following radial Hamil- 
tonian: 
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where a and J are free parameters. The potential in the 
above Hamiltonian exhibits predissociating resonances 
|16j . shown in Figjs] Because of the resonances, the ma- 
trix for this Hamihonian is non-hermitian, and so the 
time evolution operator found from U — e^*^*/'' is non- 
unitary. We will use our method described above to simu- 
late this non-unitary operator and determine the complex 
eigenvalues oi %. As a basis function for computing the 
matrix elements of the Hamiltonian, we use the orthonor- 
malized eigenfunctions of the harmonic oscillator, where 
mass equals one and the frequency equals a ref.|16|: 

vl/„(a : r) = (^f/^ -L=e''^-^' H,,{^ar). (17) 
TT v2"n! 

Hn are the Hermite polynomials of order n. For illus- 
tration purposes, we use only two bases with n — and 
n = 1 to determine the elements of the 2x2 Hamiltonian 
matrix. The matrix is found as (17': 

/1.4216 - 0.1576i 0.2782 -I- 0.2802A , , 
1^0.2782 + 0.2802i 0.6807 - 0.2361i ) ' ^ ' 

which has the eigenvector: 

, , , /-0.3790- 0.1962A 

IV',>=(^ 0.9044 j' (19) 

and the ground state eigenenergy (0.6249 — 0.4139*). The 
time evolution operator for this Hamiltonian is 

_ 0.2588 + 1.1214i -0.4569- 0.1109A (20) 
~ 1^-0.4569- 0.1109i 1.0594 + 0.7394i ) ' 

with eigenvalue 1.2268 -1- 0.8849i. From the simulation 
output, as shown in Fig|6j we find the eigenvalue of U as 
1.2259 + 0.8861i, and so the eigenvalue of V. as 0.6259 - 
0.4139i, which gives an error in the order of 10~^. This 
comes from computer round-off and the limited available 
powers of U (when the power of lambda goes to zero, 
we can no longer distinguish between |0) and |1) on the 
phase qubit; the probabilities becomes the same). The 
simulation details and the decomposition of the circuit 
are shown in the Appendix \K\ 

V. CONCLUSION 
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FIG. 5: The potential energy, V{J,a,r), as a function of 
the radial distance, r. The values of a and J are set to 

0.1. 
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FIG. 6: Output of the iterative phase estimation 
algorithm. 
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We have presented a general scheme for the execution 
of the phase estimation algorithm for non-unitary matri- 
ces. Since the circuit design used for non-unitary ma- 
trices is a general programmable circuit, the circuit for 
the phase estimation algorithm is also universal: It can 
be used for any type of matrix to compute its eigenval- 
ues. We have shown that the complex eigenvalues of non- 
hermitian quantum systems can be found by using this 
method. As an example we have shown how to obtain 
resonances for a model Hamiltonian used to describe pre- 
dissociating resonances. The method is general, and can 
be used to simulate resonances of atomic and molecular 
systems. 



Appendix A: Simulation Details 

1. The decomposition of a multi controlled network 

The uniformly controlled networks such as the one in 
FigjTaJcontrolled by k qubits can be decomposed in terms 
of 2^CN0T gates and 2^" single rotation gates[18]. For 
instance, the circuit as illustrated for fc = 2 in Fig[7a| 
can be decomposed as in Fig|7b[ The angle values in 
the decomposed circuit are solutions of the system of the 
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FIG. 7: (a) A uniformly controlled multi-qubit network, 
(b) The decomposition of the network in (a) into CNOT 
and single quantum gates. The change of the bits in the 
gray code representations determines the control qubit 

for the CNOT gates. The network in the phase 
estimation algorithm includes two consecutive circuits 
such as in (b) with and Ry single gates. 



around the z-axis, 



axis: Vf = y/Fj 



V. 



and Vjf for rotations around the y- 
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To complete this network to a whole uniformly controlled 
network, we assume that the initial four other gates are 
identity. Hence, the final decomposition for 2x2 matrix, 
the decomposed circuit includes: 4 Hadamard gates, 16 
CNOT and 16 single gates (8 CNOTs and 8 Ry gates for 
the 1//; and 8 CNOTs and 8 Ry gates for the F/), and 2 
swaps. The probability changes for the chosen state and 
the phase qubit in the simulation for this decomposed 
circuit are shown in Fig|6] 
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